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We develop a numerical scheme to make a high-frequency skymap of gravitational-wave back- 
grounds (GWBs) observed via space-based interferometer. Based on the cross-correlation technique, 
the intensity distribution of anisotropic GWB can be directly reconstructed from the time-ordered 
data of cross-correlation signals, with full knowledge of detector's antenna pattern functions. We 
demonstrate how the planned space interferometer, LISA, can make a skymap of GWB for a spe- 
cific example of anisotropic signals. At the frequency higher than the characteristic frequency 
/» = 1/(27!" L), where L is the arm-length of the detector, the reconstructed skymap free from the 
\Q • instrumental noise potentially reaches the angular resolution up to the multipoles £ ~ 10. The 

presence of instrumental noises degrades the angular resolution. The resultant skymap has angular 
resolution with multipoles £ < 6 ~ 7 for the anisotropic signals with signal-to-noise ratio S/N> 5. 

(N . 

,— I 1 PACS numbers: 04.30.-w, 04.80.Nn, 95.55.Ym, 95.75.-z, 95.30.Sf 



o 

(N 



X 



I. INTRODUCTION 



The incoherent superposition of gravitational waves coming from many unresolved sources and/or diffuse-sources 
qq ' may form a stochastic background of gravitational waves. While such a tiny signal is in nature random, their statistical 
properties contain valuable cosmological information about the astrophysical phenomena and the cosmic expansion 
f***. ■ history. Potentially, the extremely early stage of the universe can be directly probed by using the gravitational-wave 
' background (GWB), beyond the last scattering surface of cosmological microwave background. 

Aiming at detecting the tiny fluctuations, several future missions of space-based interferometers have been proposed. 
Among these, Laser Interferometer Space Antenna (LISA) will be a first-generation space-based detector launched 
in the next decade 1 . The main target of LISA is the low-frequency gravitational waves of astrophysical origin 
O" 1 around 1 ~ lOmHz. It is expected that large population of Galactic binaries produces a strong signal of GWB at 
5h ■ / lmHz Further, the extragalactic signals of GWB would be detected at high-frequency band 0,0. On the 

other hand, second-generation space interferometers, the DECI-hertz interferometer Gravitational wave Observatory 
(DECIGO) [| and the Big-Bang Observer (BBO) 2 0], have been proposed as follow-on mission of LISA. The main 
target of these missions is the primordial GWB produced during the inflationary epoch (e.g., [3, 111 llfll lll|)- For 
this purpose, the observational window around 0.1 ~ 1Hz is considered to be suitable for direct detection. With 
future technology accessible to the next two decades, direct detection will be possible for the inflationary GWB with 
amplitude ft gw ~ 10" 16 [Il[l![l|. 

The detection of GWB in the next, few decades will open a new subject of cosmology and thereby the infrastructure 
such as a new data analysis technique is needed to be exploited. Among various interesting problems, making a 
skymap of GWB is a key piece in observational cosmology. In our previous study, we have investigated the directional 
sensitivity of space interferometer to the anisotropy of GWB (Ref . . hereafter paper I). It turned out that the 
angular sensitivity of gravitational-wave detector to the anisotropic GWBs sensitively depends on the geometric 
configuration of space interferometer as well as the signal processing. Further, based on the correlation analysis, a 
method for direct reconstruction of GWB skymap has been exploited (Ref. jig], hereafter paper II). Employing the 
perturbative technique, a low-frequency skymap of GWB can be directly reconstructed (for the reconstruction in the 
low- frequency limit, see Ref. 01) • 

The aim of the present paper is to extend our previous study to the direct reconstruction of a high-frequency 
skymap beyond the low-frequency approximation. We give a simple reconstruction method, in which the problem 
is reduced to find a solution of the linear algebraic system. While the governing linear equation often becomes 
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under-determined, with a full numerical treatment, the present scheme gives an approximate but reliable solution 
for skymap reconstruction. As a demonstration, we consider a specific example for the reconstruction problem using 
the space interferometer LISA. We show that at the high-frequency band higher than the characteristic frequency 
/* = l/(27r L) ~ 9.52mHz, the LISA can make a skymap of GWB with angular resolution up to the multipoles i ~ 10. 
The presence of instrumental noises significantly degrades the angular resolution and thereby the reconstructed result 
of the multipole coefficients for GWB skymap includes a larger error. However, with the signal-to-noise ratio S/N> 5, 
the resultant skymap has angular resolution with multipoles t < 6 ~ 7. 

The paper is organized as follows. In Sec.[Hj based on the cross-correlation technique to detect the stochastic signals, 
we briefly review how to extract the information of anisotropies of the GWBs and present the simple reconstruction 
scheme, which is capable of applying to the high-frequency skymap. Several distinctions between the high-frequency 
and the low- frequency cases are discussed. In Sec. lIIII the reconstruction scheme is demonstrated in the specific model 
of the anisotropic GWB, i.e., Galactic white dwarf binaries. The quality of the reconstructed skymap is quantified in 
both the noise-free and the noisy cases. Finally, Sec. HVl is devoted to the summary and the discussion. 



II. METHODOLOGY 
A. Map-making problem 

Before addressing the reconstruction scheme, we briefly summarize how one can detect the anisotropies of GWB 
via space interferometer. 

The planned space interferometer, LISA and also the next generation detectors DECIGO/BBO constitute several 
spacecrafts, each of which exchanges laser beams with the others. Combining these laser pulses, it is possible to 
synthesize the various output streams which are sensitive (or insensitive) to the gravitational- wave signal. In particular, 
technique to synthesize data streams cance ling the laser frequency noise is known as time-delay intcrferomctry (TDI) , 
which is crucial for LISA mission ffisl IT^. |20|. see [2l| for a review). In the present paper, signals produced by TDI 
technique are used to demonstrate the skymap reconstruction. 

The multiple data streams constructed from specific combinations of laser pulses can be used to detect the GWBs 
through the correlation analysis. Suppose that one obtains the two output streams denoted by s/(t) and sj(t), the 
correlation analysis can be made by calculating the ensemble average Cu = (s i (t) s j {£)) . Naively, one may think that 
the correlation signal Cjj becomes time-independent if both the signal and the instrumental noises obeys stationary 
random process. However, this is not entirely correct, because the sensitivity of gravitational-wave detector has specific 
angular response to the GWB. Further, the motion of space interferometer is inherently non-stationary. For LISA, 
the constellation of three space crafts orbits around the Sun with a period of one sidereal year and the orientation 
of the gravitational-wave detector gradually changes relative to the sky. Hence, in presence of the anisotropic GWB, 
the amplitude of correlation signal cannot be constant, but varies in time. This is key ingredient for reconstruction 
of the GWB skymap. The basic equation characterizing the non-stationarity of the correlation signal can be written 
in the Fourier domain as (paper I, II, pi El I2II El Eg ) : 

Cij(f,t) = J ^ S h (f,n) Fu(f,Sl;t), (1) 

where the correlation signal C(f,t) is related with the one in the time-domain through the expression, Cu — 
J dfCu(f,t). The quantity Sh represents the power spectral density of GWB, which is, in general, the unknown 
function of the frequency and the sky position. On the other hand, the function Tu is called the antenna pattern 
function, which characterizes the angular response of gravitational-wave detector, and the precise functional form of 
it is known from the configuration and the orientation of the detector (see Appendix 0) . 

Eq. implies that the luminosity distribution of GWB S/,(/, fi) is reconstructed by deconvolving the all-sky 
integral of antenna pattern function from the time-series data C7j(/, t). To see this more explicitly, we decompose the 
antenna pattern function and the luminosity distribution of GWB into spherical harmonics in a sky-fixed coordinate 
system: 

s h (|/|, n) = faW/)]* ?M *) = E *tMt) Y tm (n). (2) 

£,m £,m 

Note that the properties of spherical harmonics yield p* lm — {—l) m p£^ m and a| m = (— l) l ~ m at- m , where the latter 
property comes from !Fjj(f,^l;i) = !Fij(f, —Cl;t) (paper I). Substituting into JTJ becomes 

CM*. /) = ^ E a ^(/< *) btmifW ■ ( 3 ) 

tm 
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Thus, the problem to reconstruct a GWB skymap reduces to the linear algebraic problem. That is, collecting the 
correlation signals measured at different times and solving the couples of linear equations, the multipole coefficient of 
GWB P£ m (f) can be obtained under a full knowledge of time-dependent coefficients ag m (f, t). 

Several important remarks should be mentioned in order. First, the accessible multipole coefficients pe m are severely 
restricted by the angular sensitivity of antenna pattern functions. According to paper I, the space interferometer LISA 
is typically sensitive to the multipole coefficients £ < 5 of the anisotropic GWBs in the low-frequency regime and to the 
multipoles £ < 10 in the high-frequency regime. Second, the above linear system is generally either over-constrained 
or under-determined. In the presence of the instrumental noises, this deconvolution problem tends to become under- 
determined system (paper II). As a consequence, exact reconstruction is no longer than possible in practice and we 
need to exploit an approximate method to reconstruct the GWB skymap with a limited angular resolution. 

B. Reconstruction scheme 

In what follows, assuming a prior knowledge of the time-dependent antenna pattern functions, we present simple 
reconstruction method for a GWB skymap from time-modulated correlation signals. In paper II, employing the 
perturbative expansion of the antenna pattern functions, reconstruction method for a low-frequency skymap has been 
presented. In the present paper, we especially focus on the high-frequency skymap. Here, the high-frequency means 
that the wavelength of the gravitational waves is comparable or shorter than the arm-length of the detector (or 
separation between the two detectors that produce the correlation signals), where the low- frequency approximation 
of the antenna pattern breaks down. 

The basic strategy to make a high-frequency skymap is almost the same as in paper II. Suppose that for a given 
frequency /, one obtains the discrete time-series data for correlation signals as C7j(/, U) (i — 1, 2, • • • , N). We then 
write Eq. in the matrix form as: 

c(/)=A(/)-p(/), (4) 

where the vector c has N columns, each of which contain the correlation signal C/j(/; ti). The vector p(/) represents 
the unknown multipole coefficients of the GWB spectrum and in each column, we have p\ m - Thus, if one truncates 
the spherical harmonic expansion with the multipole £ ma x, the total number of the elements in vector p(/) becomes 
(•^max + I) 2 - On the other hand, the matrix A contains the multipole coefficients of antenna pattern functions and the 
matrix has ai m (f,U) in each element. With the truncation multipole ^ max , the quantity A forms a (£ m ax + l) 2 x N 
matrix. 

As mentioned in previous subsection, the linear system Q tends to become an under-determined system, i.e., 
(f ma j + l) 2 > N, and A is generally a rectangular matrix. In such a situation, unique and exact solution of the linear 
equations cannot be obtained. In paper II, approximate treatment to solve the equations Q has been presented 
based on the idea of least-squares method. In linear system, this approximation is expressed in the following form: 

Papprox(/) = A+(/>c(/), (5) 

where A + is the pseudo-inverse matrix of Moore-Penrose type. The explicit expression is determined from the 
singular- value decomposition of the matrix A, i.e., A = W • diagK] • V. Then we have (e.g., H3) 

A+ = V f ■ diag^r 1 ] . u. (6) 

In paper II, owing to the perturbative expansion, the multipole coefficients of antenna pattern function can be 
computed analytically up to the multipole £ = 5 and the matrix A is constructed in an analytic manner. It is 
shown that the least-squares solution JSJ provides an accurate approximation for multipole moments pi m in both 
overdetermined and under-determined cases. 

In the high-frequency regime which we are interested in, the angular resolution of antenna pattern function can 
be improved and one expects that detectable multipole coefficients increase compared to those of the low-frequency 
skymap. As a price, it is no longer possible to compute the multipole coefficients ai m analytically and the the spherical 
harmonic expansion of antenna pattern function becomes a fully numerical task. In the present paper, we use the 
Fortran package, SPHEREPACK 3.1 |2^, to compute the time-dependent coefficients ai m . With a full numerical 
treatment to evaluate Eq. JSJ, the reliability of the methodology will be tested in a simple reconstruction problem 
(Sec. [HQ. 

In principle, the methodology presented in paper II can work well even in the map-making problem of the high- 
frequency GWBs. Nonetheless, to further reduce the numerical error, the least-squares method (jSJ) may be applied 
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by imposing the reality condition of the GWB intensity map, i.e., pi m (f) = {^) m P*i _™(/)- To do this, we introduce 
real quantities qi rn and rg m and divide the multipole coefficient pi m into the real and the imaginary parts as 

Pirn — qtrn + i Tim j 
Pl-m = {-I)™ {qirn ~ irim} 

for m > and 



6,0 — Qio- 



Then, Eq.© is rewritten with 



1=0 



am Ho 



+ X {aim + (-1)'™ a>l,-m} qim - i ^ {aim - (-1)™ a£,-m} ffo 



m — 1 



m— 1 



^E€;i)'A(/), 



(7) 



1=0 



where the vectors ai and pg have the (2£ + f ) columns, whose explicit expressions are 

/ a l0 \ ( q eo \ 



+ (-l) r 



7Jr 



rir, 



(8) 



where m runs from 1 to £. Hence, truncating the summation of £ by ^ max and dividing the time-ordered signals into 
the TV subsections, the matrix A and the vector p can be written explicitly 



A = — 

47T 



/ Oo(*i) ••■ ai{h) ■■■ a£ max (ti) 



\a (t N ) ■■■ ai(t N ) ■■■ ai m!lx (t N ) 



Po 



(9) 



For a reliable numerical calculation, we further divide the matrix A in Eq. @ into the real and the imaginary parts 
and apply the least-squares approximation JSJ) to the linear system ipjl. 



III. DEMONSTRATION 



In this section, general reconstruction method presented in previous section is demonstrated in the case of LISA 
detector and the map-making capability is examined for a specific example of anisotropic GWB source. 

It is expected that the low-frequency band of LISA have a strong anisotropic signals of GWB by the Galactic 
population of unresolved binaries. On the other hand, GWBs of high-frequency regime could be dominated by the 
extragalactic origin, which mainly comes from the cosmological population of white dwarf binaries 0,13- While the 
signal produced by them would reach the detectable level of LISA sensitivity, it is uncertain whether the strong 
anisotropic component exists or not. In this paper, just for illustrative purpose, we consider the Galactic GWB as a 
high-frequency source of GWB and try to reconstruct the intensity distribution of GWB. 

Fig-fflshows the projected skymap of Galactic GWB. Here, we use the simple model of Galactic GWB described in 
paper II (see also 26] ), in which we assumed that the intensity distribution of GWB just traces the Galactic stellar 
distribution observed via infrared photometry |29l| 3 . The intensity distribution shown in Fig. ^ has disk- like structure 



3 To be precise, the Galactic stellar distribution is modeled by the fitting function given in R.ef . |2fl| . which consists of the triaxial bulge 
and the disk components. 
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FIG. 1: Full resolution skymap for the simple model of Galactic GWB (paper II, The intensity distribution Sh(f, fJ) = 

P(fi) is depicted in the Galactic coordinate system. Here, the all-sky integral of intensity distribution is normalized to unity, 
i.e., JdnP{ft) = 1. 

with a strong peak at Galactic center. Spherical harmonic analysis of GWB skymap reveals that the dominant 
contribution to the intensity distribution comes from the lower multipoles with I < 4, but the contribution from the 
higher multipoles still remains significant. Even at the multipole I ~ 15, it possesses the 10% power relative to the 
monoopole component (see Table III and Fig. 12 in paper II). In this sense, it is a good exercise to diagnose the 
map-making capability of present reconstruction scheme. 

As we mentioned in Sec. Ill Al signal processing technique of LISA is called the TDI, which produces the various 
output signals canceling the laser frequency noise by combining the time-delayed six laser pulses. Among these, 
optimal TDI signals referred to as the A, E and T variables are especially suitable for the detection of GWB via 
the correlation analysis, because these outputs are basically free from the noise correlation [23, El • Here, we use the 
optimal TDIs to produce the correlation signals. The response functions for these variables are presented in the equal 
armlength case in Appendix For the reconstruction of a high-frequency skymap, the cross-correlation signals AE, 
AT and ET can be the most sensitive signals of GWB and the reconstructed skymap improves the angular resolution 
up to the multipoles I ~ 8 — 10 around the frequencies / ~ 2 — 10/*, where the characteristic frequency of LISA is 
given by /* = l/(2wL) ~ 9.52mHz. 

In what follows, specifically focusing on the frequency / = 3/* ~ 28.6mHz, we present the reconstruction results 
separately in idealistically noise-free case (Sec. IIII All and in realistic case taking account of the instrumental noises 
(Sec. ILUBl . 

A. Noise- free case 

Fig. [21 shows the intensity distribution of reconstructed skymap in the Galactic coordinate system. To obtain the 
skymap, we first create the annual modulation data of correlation signals C J 4s(/;i), C^x(/;*) and Cet(/;0 by 
convolving the original skymap Sh(f,Cl) = P (O) with antenna pattern functions in the ecliptic coordinate system. 
Here the all-sky integral of P(fl) is normalized to unity, i.e., J dtlP(tl) = 1. Dividing the one-year data of correlation 
signals into the 32 sections (N = 32), the vector c is then constructed. Combining it with the pseudo- inverse matrix 
of A, we evaluate the expression l|5|l. The resultant least-squares solution of pe m is given in the ecliptic frame and we 
finally transform it into the skymap in the Galactic coordinate system. Note that in step calculating the matrix A, 
the SPHEREPACK 3.1 package was used to compute the spherical harmonic coefficients of time-dependent antenna 
pattern functions and the multipoles higher than I = 17 were discarded (i.e., £ max = 16). 

The reconstruction result in Fig. [21 shows that due to the incomplete reconstruction, the total intensity diminishes 
and the distribution is coarse-grained, together with fake pattern. Nevertheless, the disk-like structure of the intensity 
distribution is clearly seen. Compared to the result in the low- frequency band (paper II), the angular resolution is 
greatly improved. 

It should be noticed that the present reconstruction result crucially depends on the additional parameter, i.e., the 
cutoff of singular values of the matrix A, which we call u> cu t- For a given w cu t, the pseudo- inverse matrix A + was 
constructed based on the expression JHJ using the singular- values larger than w cu t- In Fig. [2J we specifically set the 
cutoff parameter w cut to 10 -7 . The cutoff value directly affects the resolution of reconstructed skymap. To show the 
significance of this, in Fig. |3 reconstruction results of skymap with various cutoff are plotted. Clearly, larger cutoff 
value makes the angular resolution of skymap worse. It seems favorable to use the small value of u> cu t- However, as 
indicated in the expression JfjJ), the result using the small singular- value tends to be sensitively affected by the errors 
associated with numerical computation of pseudo-inverse matrix and/or the instrumental noises. Hence, we need to 




FIG. 3: Dependence of the reconstructed results on the cutoff parameter of singular values in absence of instrumental noises: 
uicut = I0~ e (left), 1CT 4 (middle) and 1Q~ 2 (right) . 



seek an optimal value of w cu t to get a better reconstruction result. This will be discussed in next subsection. 

Now, we wish to quantify the quality of the reconstructed skymap and discuss the validity of present reconstruction 
scheme. To do this, we introduce the correlation parameter defined by 



^*corr (-^cut ) 
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In the expression Hl()|) , while the quantity 5 , ^ reconst ) represents the reconstructed skymap from the least-squares solution 
(|5]l , the function 5^ true ) j s the skymap of the true intensity distribution dropping the higher multipole moments with 
& > lout- To be precise, this is defined by 
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The correlation parameter r COTr (£ cut ) quantifies the degree of similarity between the reconstructed skymap and the true 
skymap subtracting the higher multipoles. Hence, evaluating r corr (^ cut ), one can characterize the angular resolution 
of the GWB skymap. In addition to this, we also introduce the averaged fractional error for the multipole coefficients 
Err [pi 

rn ], given by 



Err [p e „ 
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(true) 
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(13) 



which quantifies the accuracy of the determination of multipole coefficients. With the two quantities f C orr(^cut) and 
E rr [pim], the quality and the accuracy of the reconstructed skymap can be quantified. 

Left panel of Fig. 0] shows the correlation parameters for the reconstructed skymap with various cutoff values 
w C ut- For a fixed value ui cu t, the correlation parameter as function of truncation multipole £ cu t has a single peak 
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and the peak value is typically r corr ~ 0.92. As decreasing the cutoff value u> cu t, the location of peak tends to be 
shifted to a larger £ cu t, indicating that the angular resolution becomes improved. For small cutoff value w cu t = 10 -6 , 
reconstructed skymap is qualitatively similar to the true skymap dropping the higher multipoles £ > 12. 

Right panel of Fig. 0] shows the root-mean-square value of the averaged fractional error of pi m , he., Err[p^ m ]. 
As anticipated from left panel, the smaller value of w cu t tends to decrease the errors in each multipole coefficient. 
However, due to the low-resolution of the antenna pattern function, the reconstructed results for higher multipoles 
I > 10 become almost vanishing. As a result, the fractional error Err [pi m ] of £ > 10 approaches unity. It should 
be noticed that the fractional error also becomes unity for the monopole moment of GWB signal. This reflects the 
important properties that the cross-correlation of optimal TDI variables are generally blind to the monopole intensity 
for the AT and ET signals and also to the dipole anisotropy for the AE signal. It has been shown in paper I that 
these features generally hold irrespective of the observed frequency band. In this respect, map-making issue with 
LISA would be problematic in determining of the monopole intensity. Apart from this, right panel of Fig. 0] implies 
that the present reconstruction scheme has a potential to provide a precise determination of multipoles. With the 
cutoff value w cu t — 10~ 6 , the accuracy can reach less than a few percent level. Of course, these are the outcome based 
on the idealistic situation free from the instrumental noises. In next subsection, the significance of the instrumental 
noises on the map-making problem will be clarified. 



B. Noisy case 



In reality, the instrumental noises are additively mixed into the time-series data of gravitational-wave signals. While 
the cross-correlation variables AE, AT and ET themselves are free from the noise correlation, the correlation analysis 
with finite number of samples is inherently affected by the statistical fluctuations, among which the instrumental 
noises becomes the most dominant component in the weak-signal case. To discuss their influences on the map-making 
problem, we first quantify the signal-to-noise ratio for the stochastic GWB. Based on the cross-correlation statistic 
with a suitable filter function, the signal-to-noise (S/N) ratio for the stochastic GWB in the specific frequency interval 
/ - A//2 ~ / + A//2 is defined by (e.g., MM) 



where the summation (I, J) represents the sum over the whole cross-correlation data (i.e., AE, AT and ET). The 
quantities T or bit and T b s respectively denote the orbital period of LISA corresponding to the one year and the 
total observation time longer than the orbital period. Below, we specifically set the parameters to T b s = 10 8 sec and 
A/ = 10~ 3 Hz. The function Nj(f) represents the noise spectral density for the output data /. The explicit functional 
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FIG. 5: Annual modulation data of cross-correlation signals Cae (left), Cat (middle) and C'et (right) measured at the 
frequency / = 3/* ~ 28.6mHz. Top and bottom panels show the results with signal-to-noise ratio S/N= 5 and 30, respectively. 
In each panel, filled and open circles represent the real and the imaginary parts of the noisy signals Cu generated according 
to Eq. 1151 . While the continuous thin lines indicates the correlation signals free from the detector noise, thick lines indicate 
the fitted result of the noisy signals to the harmonic functions llbl . 



form of the noise spectrum is given in Appendix ^ together with the specific parameters of the proof mass and the 
optical-path noises for LISA. 

We are specifically concerned with the qualitative change of the reconstructed skymap in the presence of instrumental 
noises. For this purpose, rather than performing a large-scale extensive simulation that mimics a realistic signal 
processing [33l l34l l35| , we here perform a very simple simulation in which the random noises are added by hand to 
the (noise-free) correlation signals computed in previous subsection. To be precise, for z-th sub-sectional data of the 
totally N = 32 correlation signals Cu, we generate the time-series random data Cu as 

c ^m = c I Af;k) + [ {Toba/N)Af j (* = !,•••, AO (15) 

where the variables are the Gaussian random variables with zero mean and unit variance, i.e., = and 
= 1. In the above expression, the first term in right-hand-side is the theoretical time-modulation signals used in 
previous subsection, where the amplitude of GWB signal, Sh(f, fi) = A ■ P(£l). On the other hand, the second term 
represents the randomness arising from the instrumental noises, whose amplitude is estimated based on the S/N ratio 
in weak-signal limit. 

Fig. shows the time-series data of cross-correlation signals for Galactic GWB in the case of S/N= 5 (top) and 
S/N= 30 (bottom). The open and the filled circles represent the real and the imaginary parts of the noisy data Cu, 
respectively. Note that in plotting the data, we set A = 1 so that the all-sky integral of GWB spectrum is normalized 
to unity and the amplitude of the random noises is appropriately rescaled according to the S/N values. As it is clear, 
direct use of the raw noisy data makes the quality of the reconstructed image worse. Hence, we tried to fit the noisy 



FIG. 6: Reconstructed skymap from the noisy data (see Fig. 0. Left and right panels respectively show the results from the 
S/N=5 and 30 noisy data. From top to bottom, the cutoff parameter of singular values were set to iw C ut = 10 -3 , 5 x 10~ 3 , 
1(T 2 and 5 x 10~ 2 . 

data Cu to the harmonic functions fu(t): 

fcmax 

fu(t)= c k e lk ^ 1 - Worbit = 27r/T orbit . (16) 

fc=-fe max 

The resultant fitting functions were then used to perform the reconstruction scheme. In Fig. [5] the thick-solid and 
the thick-dashed lines are the results fitted to the harmonic function (|16f) with fc max = 8. For comparison, we also 
plot the continuous thin lines as the noise-free correlation signals. 

Fig. H3 shows the reconstructed results from the noisy data for the various cutoff parameters of the singular- values, 
Went'- Went = 10 -3 , 5 x 1CT 3 , 10 -2 and 5 x 10 -2 from top to bottom. The left panel plots the results in the S/N= 5 
cases, while the right panel shows the skymap in the S/N= 30 cases. The skymap with small cutoff value w cu t = 10~ 3 
is affected by the instrumental noises and the resultant intensity map show featureless fake patterns. As increasing 
the cutoff parameter, fake intensity pattern gradually disappears and the strong intensity peak seen in the original 
GWB map becomes prominent. The quality of the resultant skymap with large u> cut depends on the signal-to-noise 
ratio, S/N. For the cutoff parameter u> cut = 10~ 2 , the resultant skymap in the S/N=30 case is very similar to the one 
in the noise-free result with u> cut = 10 -2 (right panel of Fig. while the fake intensity image is still dominant in 
the S/N=5 case. 

The influence of instrumental noises shown in Fig. [5] may be easily deduced from the expression © with ©. In the 
presence of the noise, the vector c which represents the correlation signals additionally contains the noise term n and 
one can write it as c = A ■ p + n. Then, the least-squares approximation (J5J leads to the unwanted term A + • n. Since 
the pseudo- inverse matrix A + contains the reciprocal of the singular values, {u^ 1 }, this additional contribution can 
become dominant and affect the final reconstructed result unless introducing a larger cutoff value w CU f Therefore, in 
order to reduce the influence of the noise term, the cutoff parameter should be, at least, set to ^(S/N) -1 . 
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FIG. 7: Cross-correlation parameters of reconstructed skymap in presence of the instrumental noises with S/N=5 (left) and 
S/N=30 {right). In both panels, the cutoff parameters of the singular- values of the matrix A were set to w C ut = 10 -3 (solid), 
5 x 10 -3 (long- dashed), 10 -2 (short- dashed) and 5 x 1CP 2 (dotted), from bottom to top. The error bars are estimated from the 
sample variation among 100 realizations. For comparison, the cross-correlation parameter in the noise-free case (w cut = 10 -6 ) 
is also plotted in thick solid line. 



In Fig. \7\ quality of the reconstructed skymap is quantified by evaluating the correlation parameters r corr . Also, 
the averaged fractional errors Err [pe m ] are computed and are presented in Fig. [S] for various cutoff values. In both 
figures, the error bars indicate the la variation among 100 realizations. As anticipated from Fig. the quality of 
final skymap is significantly degraded in the case of the low S /N data. The result with a smaller cutoff value w cut has 
a little correlation with the true skymap due to many fake patterns, which is mainly attributed to the errors in the 
higher multipoles I > 8. For the high-signal-to-noise ratio S/N=30, the situation is somehow improved. With a cutoff 
value around u> cu t ~ (S/N) -1 , the quality of reconstructed image becomes similar to the true map with £ < 6 — 7. 
The fractional error in each multipole coefficient reaches ~ 10%, although it seems still miserable compared with the 
noise-free result. 

Finally, accuracy of the reconstructed multipoles achieved with a given signal-to-noise ratio S/N is shown in Fig. [5] 
Except for S/N=100, the cutoff value was all fixed to w cu t = 5 x 10~ 2 . For high signal-to-noise ratio S/N=100, the 
accuracy is further improved adopting the small cutoff parameter, u> cu t = 5 x 10~ 3 . In Fig. EI apart from the monopole 
moment and the higher multipoles of t > 8, there exist some multipoles that are still difficult to determine. Recalling 
from the similar trend found in the noise- free case (right panel of Fig.^J), the non- uniformity of the accuracies might 
be attributed to the angular response of the antenna pattern functions at the frequency / = 3/*, not the character 
of the specific model of anisotropic GWB. This implies that the multi- frequency data analysis is important for the 
accurate determination of the multipole moments pi m . Anyway, the quantitative estimation of p£ m at few percent 
level seems very difficult even when S/N=100. This means that in contrast to the low-frequency case, the present 
reconstruction method might not be so powerful to determine each multipole moment of high-frequency skymap, but 
should be rather helpful to study the all-sky distribution of anisotropic GWBs in bird-eye view. 



IV. SUMMARY & DISCUSSION 



In this paper, we have discussed the reconstruction of a high-frequency skymap of gravitational-wave backgrounds 
with a space-based interferometer. Owing to the cross-correlation technique, we present a simple numerical scheme to 
reconstruct a skymap of anisotropic stochastic signals. While the methodology presented here is basically the same 
one as described in paper II, i.e., least-squares solution of linear algebraic system, there are several distinctions in 
the reconstruction of the high-frequency skymap. First of all, the antenna pattern function becomes complicated and 
the analytic treatment based on the perturbative expansion of antenna pattern function is intractable. Hence, a full 
numerical treatment to reconstruct the higher multipole moments is necessary. On the other hand, the sensitivity to 
the high-frequency signals is improved and no multi-frequency data is needed for the reconstruction of GWB skymap. 
As a result, the skymap with LISA demonstrated in the Galactic GWB case yields a better angular resolution up 
to the multipoles I ~ 10 in the absence of the instrumental noises. The presence of noises degrades the angular 
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FIG. 9: Same as in Fig. |^1 but the dependence of signal-to-noise ratio is shown. The solid, long-dashed, short-dashed and 
dotted lines respectively show the results with S/N=100, 30, 10 and 5. In each case, except for S/N=100, the cutoff parameter 
w cu t was set to 5 x 10 -2 . In the case of the high signal-to-noise ratio S/N=100, fractional error can be further reduced if one 
adopt a smaller cutoff value. In this plot, the result with ui cu t = 5 x 10~ 3 is shown. 



resolution and the multipole coefficients for GWB skymap includes a larger error. However, with the signal-to-noise 
ratio S/N> 5, the resultant skymap has angular resolution with multipoles I < 6 ~ 7, which gives a better quality 
compared to the one achieved by the low-frequency skymap. 

Since the methodology presented here is still very primitive, there would be several extensions to improve the 
angular resolution of reconstructed skymap. Perhaps, one naive approach is the combination of the multi-frequency 
data set. With the multiple data, one effectively increases the signal-to-noise ratio. Also, the directional information 
for anisotropic signals would be obtained additionally, through the frequency-dependent angular response of antenna 
pattern functions. Another possible approach is to combine the parametric models of GWB distribution characterized 
by the finite number of parameters and to determine these parameters. With fewer parameter set, the present 
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methodology tightly constrains the model parameters. 

The reconstructed skymap obtained from the present scheme may be regarded as the first step of the synthesis 
imaging processing, like a dirty map in radio astronomy |3(| . Under a priori information about the source distribution 
and the reasonable assumptions, it would be CLEANed by using the iterative deconvolution algorithm (see Ref.[37j 
for the application of CLEAN algorithm to the data analysis of gravitational waves). That is, with the compact GWB 
distribution coming from the nearby sources, one can create the synthesized GWB map convolving with the antenna 
pattern functions. We then subtract each compact component from the dirty map and repeat the procedure iteratively 
until all significant source structure has been removed. A CLEAN map is finally obtained from the residual intensity 
distribution by adding the removed GWB components with suitably smoothed sampling function. This CLEAN map 
would be useful and helpful to discriminate the GWBs of nearby origin and cosmological origin and even to identify 
the specific GWB signals. 

Finally, the map-making problem considered in the paper is related to the tomographic reconstruction technique to 
resolve the distribution of the binaries (3^. While the present technique only relies on the amplitude of the signals, 
the tomographic approach fully takes account of the phase information. As a result, the angular resolution is greatly 
improved and the identification of signals becomes efficient even in the crowded samples. This may be a hint to 
improve the angular resolution of GWB skymap. 
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APPENDIX A: ANTENNA PATTERN FUNCTIONS AND INSTRUMENTAL NOISES 

Here, we give an explicit expressions for antenna pattern functions for LISA used in the main text. First recall the 
definition of the antenna pattern function (paper I, II, |24j): 

T u (f,Cl; t) = e^/fi.Cx.-x,) Ff*(nj; t)Ff(tl,f; t) 



A=+,x 



F/(f2,/; t)=Dj(n,/; t):e A (n), 



(Al) 



where e A is the polarization tensor for GWB and D/ is the detector tensor for each output signal, to which we 
specifically adopt the optimal TDI variables A, E and T. The advantage of using the optimal TDIs in the cross- 
correlation analysis is that the noise spectra for cross-correlation data becomes exactly vanishing [sH Isij . These 
variables are simply realized by combining the three Sagnac signals called Si, S2 and S3 (often quoted as a, (3 and 7 
in the literature): 



D A ^ 7! (D S3 



DsJ, 



De = -^(D Si -2D S2 



Dn 



+ Ds 3 ), 
Ds 3 ). 



(A2) 



The detector tensor for Sagnac signals can be obtained from the time-delayed combinations of one-way Doppler 
tracking calculations for optical-path length [lja, |3j| |4jJ. For example, the Sagnac signal Si measures the phase 
difference between two laser beams received at space craft 1, each of which travels around the LISA array in clockwise 
or counter-clockwise direction. Then, in the equal-arm length limit, the detector tensor for Si becomes 

D Sl (n,/;t) = ^[{a(t)®a(t)}T a (n,/;t) + {b(t)®b(t)}T b (n,/;t) + {c(t)®c(t)}a;(n,/;t)] ; 



T a (n,/;t) = e - 3 ^/ 2 ^e-™{- 2+a W-">sinc 



(l + a(t)-n) 



T b (fi,/; t) 



e -i(//2)[3+{a(t)-c(t)}-n] I ginc 



(i + b(t)-n) 



e - l (//2){2+a(t).0} sinc 



(l-b(t)-n) 



/ 



(1 - a(t) ■ n) 
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%(Sl,f; t) 



= e~ 3if / 2 



e -i(/72){2-c(t).n} sinc 



|(i + c(t)-n) 



e «(//2) {2+ c(t).ri} sinc 



c(t) • SI) 



where the time-dependent vectors a(t), b(t) and c(t) represent the unit vectors pointing from the space craft 1 to 2, 
2 to 3 and 3 to 1, respectively. Here, the quantity / denotes the normalized frequency defined by 



/ 



f_ 

/* 



/, = 1/(2ttL). 



(A3) 



With the arm-length of LISA L — 5 x 10 6 km, the characteristic frequency /» becomes 9.52mHz. The analytic 
expressions for other detector tensors Dg 2 and Dg 3 are also obtained by the cyclic permutation of the unit vectors a, 
b and c. 

In the expression of antenna pattern function (|A1|) and/or detector tensor, the time-dependence is incorporated 
through the unit vectors a(t), b(t) and c(t), which are determined by the orbital motion of LISA. Assuming the rigid 
motion in the equal- arm length limit, these are given by 4 




R*(-M*)) • Ry(-fe) • R*(M*)) 




(A4) 



Here the quantities R y and R 2 are respectively the 3x3 rotation matrices around y- and z-axes defined in the ecliptic 
coordinate system (see Fig.l of paper I). The angles 4>d and Qd are chosen as <f>jj = —n + cij or bit i and 8o — — tt/3. 
The vectors ao, bo and Cq represent the orientation of LISA in detector's rest frame. Here, we specifically set 



a 



V3 1 



b = (0, -1,0), c = 



V3 1 



(A5) 



Fig. llOl shows the intensity distribution of the antenna pattern functions at t = 0. The frequency is specifically chosen 
as / = 3/* ~ 28.6mHz. The antenna pattern functions at high-frequency regime give a different directional response, 
which implies that the cross-correlation signals AE, AT and ET mutually provide an independent information about 
source distribution. Note that the typical angular size of the intensity patterns is roughly 30° ~ 60°, which basically 
limits the angular resolution of the final reconstructed skymap. 

Finally, the noise spectral density for each output signal is presented, which is necessary to estimate the signal- 
to-noise ratios. Although the optimal TDI variables are free from the noise correlations, non-vanishing contribution 
to the self-correlation signals does exist. The noise spectral densities for optimal TDIs are calculated as (see also 
El El): 



N A (f) = N E (f) = sin 2 (//2) {8 (2 + cob/) S shot (f ) + 16 (3 + 2 cos / + cos 2/ ) ,s' a 
N T (f) = 2 (l + 2cos/) 2 {5 shot (/)+4sin 2 (/72) 5 accc i(/)}, 



i(/)}, 



(A6) 



where the variables S s hot and 5 acco i represent the proof mass noise and the optical path noise, respectively. Adopting 
the numerical values reported in Ref.pLTI] (see also Hjj), we obtain 5 s hot(/) = 1-60 x 10~ 41 Hz -1 and 5 acc el(/) = 
2.31 x 10- 41 (mHz//) 4 Hz- 1 . 
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